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Abstract. The aim of this paper is the accurate numerical study of the KP 
equation. In particular we are concerned with the small dispersion limit of 
this model, where no comprehensive analytical description exists so far. To 
this end we first study a similar highly oscillatory regime for asymptotically 
small solutions, which can be described via the Davey-Stewartson system. 
In a second step we investigate numerically the small dispersion limit of the 
KP model in the case of large amplitudes. Similarities and differences to the 
much better studied Korteweg-de Vries situation are discussed as well as the 
dependence of the limit on the additional transverse coordinate. 
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1. Introduction 

This work is concerned with the 2 + 1 dimensional Kadomtsev-Petviashvili equation 
(KP), given by 

(1.1) d x (d t u + ud x u + e 2 d xxx u) +\d yy u = 0, A = ±1, 

where (t, x, y) £ M. t x K x x M. y and where e > is a (small) scaling parameter, 
as introduced below. The case A = — 1 corresponds to the so-called KP-I model, 
whereas A = +1 is usually referred to as KP-II equation. Both cases can be de- 
rived as models for nonlinear dispersive waves on the surface of fluids [21] (see 
also [22] )• Thereby it is assumed that the propagation of the waves is essentially 
one-dimensional with weak transverse effects. Roughly speaking, KP-I describes 
the case when surface tension is strong, whereas KP-II is a good model for weak 
surface tension. Moreover, the KP type equations also arise as a model for sound 
waves in ferromagnetic media |42| and in the description of two-dimensional non- 
linear matter- wave pulses in Bose-Einstein condensates, see e.g. |19l I23| . 
Notice that the KP equation has originally been derived in the following form 
d t u + ud x u + e 2 d xxx u + \d y v = 0, u\ = ui(x,y), 

(1-2) 

a y u = d x v, 

where the second equation can now be seen as a constraint for the Cauchy problem 
given in the first line. Therefore we shall also impose for equation (|1.1|) initial data 
in the Schwartz class of smooth and rapidly decreasing functions, i.e. 

(1.3) u| t=Q = ui G 5(M X x %), 
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even though an interpretation of in terms of a Cauchy problem seems not 

completely obvious at first sight, due to the appearing mixed derivative d x tu. Nev- 
ertheless this is the usual approach when dealing with l|l.lfl and it is strengthened 
by the fact that in large parts of the literature one considers the KP model in its 
so-called evolutionary form 

(1.4) d t u + u d x u + e 2 d xxx u + A d~ 1 d yy u = 0, u| i=Q = u\{x, y). 

Here, the anti-derivative 9" 1 can be uniquely defined via 

(1-5) d^fix) := \ (J ^ /(C) dC - /(C) dc) , 

at least if fix) decays sufficiently fast to zero, as x — » ±oo, respectively. Alterna- 
tively ftjT 1 can be seen as a Fourier multiplier with the singular symbol —\/k x , as 
we shall do in the following. 

In this paper we are interested in the accurate numerical simulation of the KP 
equation using a spectral scheme with preconditioning. To this end we remark that 
alternative numerical schemes for 2 + 1 dimensional wave equations have been pro- 
posed in [111 \V2\ . Moreover, earlier numerical studies of the KP equation can also 
be found in, e.g., 1210 EH E] QH] , aiming mostly at the description of (interact- 
ing) solitons. In our work though the focus will be mainly on asymptotic regimes 
of the KP model. In particular we are concerned with the so-called dispersionless 
limit, i.e. the limit of l|l.l|) as e — > 0. 

To motivate this study consider the unsealed (dimensionless) KP model 

(1.6) dx(dru + udxu + dxxxu) + Xdyyu — 0, u\ T=0 = ui(X,Y), 

where (T, X, Y) are now considered as the natural scales of observation, or micro- 
scopic scales. In ljl.61) we shall now introduce slowly varying solutions of the form, 
i.e. 

(1.7) u = u(eT,eX,eY), < e < 1, 

where e is a small scaling parameter, i.e. the microscopic/macroscopic scales ratio. 
These slow variations, for example, can be introduced via a class of corresponding 
initial data. Plugging (|1 .T|) into (|1.6J) and denoting (t,x,y) — (eT,eX,eY), we 
obtain the scaled equation (f 1 . 1|> with e <C 1. In the following we shall always 
consider the KP equation in the form (f 1 . 1 1> . The notion of the coordinates thus 
refers to the macroscopic sales. 

In the formal limit e->0we get the dispersionless KP equation (dKP), i.e. 

(1.8) d x (d t u + ud x u) + XdyyU = 0. 

Clearly this is a singular limiting procedure, which requires particular care. The 
analogous problem for the scaled Korteweg-de Vries equation (KdV), see equation 
(|2.1(l below, has been intensively studied by several authors, cf. [131 12"§1 13T| . 
More recently this problem has also been treated numerically in In the KdV 
case, it is known that the solution of the Cauchy problem for smooth hump-like ini- 
tial data is characterized by the appearance of a zone of rapid modulated oscillations 
|16H39| . A rigorous asymptotic treatment of these oscillations has been established, 
relying heavily on inverse scattering techniques and complete integrability. On the 
other hand, a complete mathematical description of the small dispersion limit for 
the KP model has not yet been achieved EHj • The purpose of this paper is 
thus to explore this problem numerically. Similar to KdV, it is expected that the 
dKP model will develop shocks in finite time. In this case, the formal limit (|1.6|) is 
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no longer an accurate description of as e — > 0. For small, but still non-zero 

e, the dispersive term oc d X xxU in ljl.4|) presumably will smooth out these shocks 
by forming an oscillatory zone. This work provides numerically evidence for these 
facts. 

More precisely the outline of this paper is as follows: In the next section we shall 
first recall some well known features of the KP equation. Also, we shall derive 
the Green's function of the corresponding linear model, which already gives some 
insight on the oscillatory nature of the considered equation. In Section the nu- 
merical algorithm is presented and we shall test its performance for two cases of 
soliton propagation. In Section0]we discuss the the asymptotic behavior for e<l 
in the particular case of small solutions, i.e. with amplitudes of order 0(e). This 
yields an asymptotic description via the Davey-Stewartson system, which is studied 
numerically for the case A = 1. The corresponding asymptotic error to the KP-II 
equation, defined in the L 2 (M. 2 ) sense, is, within numerical precision, found to be 
0(e 5 / 2 ) which is remarkable since it is exactly as expected from the corresponding 
rigorous KdV studies in |35j . Finally we turn to the problem of small dispersion 
for solutions of order 0(1) in Sectional We shall numerically investigate the en- 
tropy solution of the dKP equation and the corresponding oscillatory behavior in 
the KP model, in particular its dependence on the y-coordinate. We also find nu- 
merically that the asymptotic L 2 (M 2 )-error in this case is 0(e 3 / 2 ), at least before 
the appearance of the first shock in the corresponding dKP solution. 



2. Properties of the KP equation 

2.1. Preliminaries. In the context of water waves, the KP equation can be seen 
as a two-dimensional generalization of the celebrated KdV equation i-e. 

(2.1) d t u + ud x u + e 2 d xxx u = 0, u\ t=Q =Ui(x). 

Clearly, any solution to KdV is a y-independent solution to the KP equation, form- 
ing the KdV sector of (jl.l|) . Indeed similar to KdV, the KP model belongs to the 
class of integrable systems and is known to conserve an infinite number of quantities, 
among them mass 

(2.2) M[u{t)]:= u 2 (t,x,y)dxdy 

JR 2 

and energy 
(2.3) 

E[«(t)] := ^d x u(t, x. y)) 2 - X(d~ 1 d y u{t, x, y)) 2 - j u 3 (t, x, yfj dxdy, 

as well as momentum etc., cf. pQ for more details. Again similar to KdV, explicit 
solutions to the KP model can be obtained via inverse scattering methods. In the 
following though we shall not go into more details on this important topic, arising in 
the study of integrable systems, and rather refer to, e.g., OESj and the references 
given therein. 

Despite the apparent similarity of KP-I and KP-II, the two equations behave rather 
differently both qualitatively and quantitatively. For example, only the KP-I has 
the advantage of having an energy which is positive definite in the leading order 
terms. Additionally, unlike the KP-I model, the KP-II equation does not admit 
so-called lump solitons. Generally speaking, a KP soliton is a solution of the form 



(2.4) 



u(t, x, y) = u(x -tv x ,y- tv y ), 
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with (constant) group velocity v g — (v x ,v y ) G R 2 . Clearly, any KdV soliton is also 
a KP soliton solution. Wc also remark that both KP models bear several differences 
in what concerns the stability or instability of solitons, cf. [SJI^ and the references 
given therein. The single lump soliton to ijl.lfl , with A = — 1 , has the following form 
(here, we set e = 1 for simplicity) 

n 2Ac(l-c(x-3ct) 2 + 3c 2 y 2 ) 

(2 - 5) (iT^WWr ' ceM - 

Note that the lump soliton only decays algebraically in both spatial directions, 
as \x\, \y\ — > oo. On the other hand both equations admit so-called line soliton 
solutions. These are solitons which are infinitely extended in one spatial direction 
and which exponentially decay in the other direction. The simplest solution of this 
kind is a y- independent 1-soliton of the KdV equation (again e = 1 for simplicity) 

12 

(2.6) u(t, x) = § i with xq G K fix. 

cosh (x — xo — At) 

Solutions which are even better localized in x and y cannot be expected (see also 
[7]). Indeed it is a common feature of both KP-I and KP-II that, even if the initial 
data is a Schwartz function the corresponding solutions u(t) for t > in general 
will not stay in the Schwartz class, unless u\ satisfies in addition an infinite number 
of constraints. Rather it is known that, for generic initial data, the solutions u(t) 
will develop "tails" decaying only algebraically in certain spatial directions. The 
reason for this behavior is already seen on the level of the Green's function for the 
linear KP equation, as we shall discuss in the next subsection. 

2.2. The linear KP equation. In the following we shall study the linear part of 
the KP equation, i.e. 

(2.7) d x (d t u + e 2 d xxx u) + Xdyyu = 0, u| t=Q = ui(x,y). 

for smooth (finite mass) solutions u(t) S L 2 (M. 2 ) n L 1 (M 2 ). If we denote by 

{Tu{t)){k x ,k y ) = u{t,k x ,k y ):= f u(t,x,y)e^ x+k yy'>dxdy, 



the Fourier transform of u(t) w.r.t x and y, the above given equation (|2.7|) is 
equivalent to 

(2.8) ik x d t u + e 2 k A x u-Xk 2 y u = 0. 

In order to pass to the evolution form (|1.4f> we need to be able to apply the singular 
multiplier —\jk x to this equation, to obtain 

( \k 2 \ 

(2.9) d t u + il-^-e 2 k 3 x ju = 0. 

Having in mind that u(t) £ Cb(K. 2 ), for u(t) G i 1 (M 2 ), this certainly requires 
that u(t,0,k y ) = 0, for all k y G M. To ensure this, we shall impose the following 
condition on the initial data 

(2.10) f Ul (x,y)dx = 0, 

which, in Fourier space, is equivalent to ui(0, k y ) — 0, for all k y G M. In this case 
the (weak) solution of (|2.9|) is obviously given by 

(2.11) u(t, k x ,k y ) = e- it{xk l/ k *-^ Ul {k x ,k y ), 
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which satisfies u(t,0, k y ) = 0, for all t £ R. The constraint 1(2.10(1 thus seems to be 
quite natural when passing to the formulation 1(1.4(1 , cf. EH] ■ Moreover it can 
also be seen from the numerical example given in Section |3 that initial data which 
do not satisfy 1(2.10(1 become non-differentiable for any t ^ 0. 

Remark 2.1. Alternatively one might also remark that (formally) integrating p. 1(1 
w.r.t. to x £ M. gives 

(2.12) / d yy u(t, x, y) dx = 0, Vyel,t^0, 
Jr 

or, equivalently, — fc 2 u(t, 0, k y ) = 0, a condition certainly fulfilled by l|2.11|l . 1(2.10(1 . 

Of course u(t,0,k y ) = is not sufficient to obtain a strong KP-solution (being 
differentiable in-time) which satisfies 1(2. 9|l point wise in Fourier space. To this end 
one would require a sufficient fast decay of u(t, k x ,k y ), as k x — > 0. Therefore most of 
the literature is concerned with the mild formulation ((2.11|) . which is consequently 
translated to the non-linear model, via Duhamel's formula 

(2.13) u(t,x,y) = U(t)ui(x,y) - U(t - s)u(s,x,y)d x u(s,x,y)ds. 

Jo 

Here, U(t) denotes the unitary group in L 2 (M. 2 ) defined via its symbol in Fourier 
space, i.e. 

(2.14) U(t) := exp (-it{\k 2 y /k x - e 2 k 3 x )) . 

Following this approach, well posedness issues are usually studied in non-isotropic 
Sobolev spaces H Sl (R x ) x H S2 (R y ). It should be noted that local and global well 
posedness results are much more complete in the KP-II case than in the case of 
KP-I. The first result in this direction has been obtained in [8] proving, that the 
KP-II model is globally well posed for u\ £ L 2 (M 2 ), or u\ £ L 2 (T 2 ). As far as 
Sobolev regularity is concerned we shall only remark that the required index pair 
for KP-II is: si > -1/3, s 2 > [S3, referring to [El E3 EH1 ESj for more details on 
these issues. 

Remark 2.2. Note that the formal identification d tx u = d xt u requires particular 
care if one does not consider solutions, which are differentiable in time |32| . In this 
context we also remark that sometimes even weaker notions of solutions which are 
only defined to be distributional in-time, are taken into account, cf. UJ. 

In order to get more insight on the dynamics of the linear equation we first note 
that 1(2.11(1 implies 

(2.15) u(t, x, y) = ^ (j^Cty)) * ui(x, y), 

where * denotes the convolution w.r.t. 16I and y £ ML In order to calculate the 
inverse Fourier transform of U (t) we shall split the transformation in two parts via 

(2.16) T~ x L-^KI**-^S\ = T~ x ( e - itAfc ^) * T- x (e"^) . 

In the following let us restrict to the case of solutions for times t > 0, for simplic- 
ity. Next, consider the last factor on the r.h.s. of ((2.16(1 and recall the Fourier 
representation of the so-called Airy function, i.e. 

(2.17) Ai(x) = — — f c i{xk+k3/3) dk. 
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It is well known, see, e.g., 0, that the function Ai(x) admits a series expansion in 
the form 

Moreover Ai(x) obeys two completely different asymptotic limits as x — > ±00, 
respectively, namely 

. , X-.+00 l -^11 1 /2a; 3 / 2 tt 

(2.19) Ax(x) ~ e 3 ' M*) ~ ^T7l cos ^— - ^ 

We consequently obtain, by substituting k — k x {ite 2 ) 1 /' i , that 

(2.20) ^ (e«^) = ^ Ai (^) ® v* > 0, 

in S'(M. X x M. y ), i.e. the sense of tempered distributions. 

On the other hand, computing the inverse Fourier transform of the first factor on 
the r.h.s. of i|2.16() yields, after some lengthy computations, that (for t > 0) 

!4£l/2 
else. 

Again this has to be interpreted as a distribution on Schwartz space. (We note 
that formula (|2.21|) is equivalent to the one found in [Sj, where only the KP-I case 
is considered though.) Expression 12.2111 yields the Green's function for the linear 
partial differential operator 

(2.22) Lf := (d xt + Xd yy ) f = 0, A = ±1. 

Note that Lf = is a hyperbolic PDE as can be seen by introducing the coordinates 
a — x + t and (3 — x — t, and transforming l|2.22|l into 

(2.23) (d aa ~dpp + \d vv )f = 0, A = ±l. 

This is nothing but the standard 2 + 1 dimensional wave equation, where a or (3 
take on the role of a "time variable" for A = —1 and A = 1, respectively. Since 
the operator L can also be interpreted as the linear part of the dKP equation 
()1.8(l , imposing at time t = an initial condition of the form i|1.3fl indeed furnishes 
a characteristic Cauchy problem. Moreover, the space-time region determined by 
Axt + Xy 2 — is also a characteristic which explains the divergences of the Green's 
function 12.21jl there. 

Combining l|2.15(l with the formulas H2.20JI and (|2.21(l we can draw several important 
consequences: 

• The y-dependence of the solutions to the linear KP equation (|2.7(l is com- 
pletely non-oscillatory and has a rather slow algebraic decay as \y\ — > 00. 

• For any time t>0we obtain infinitely extended tails of the solution within 
the region determined by Axt + Xy 2 > 0, for t > 0. 

• The oscillatory behavior of the linear KP equation is governed by the Airy 
function. In our scaling this yields oscillations with wave length 0(e) in 
the ^-direction, as can be seen from H2.19|l . 
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3. The numerical algorithm 

3.1. A Spectral approach with preconditioning. We are interested in the 
numerical solution of the KP equation for rapidly decreasing, smooth initial data U\. 
To study the zone of fast oscillations, due to e <C 1, it is convenient to use a periodic 
setting of sufficiently large period. This allows for the use of spectral methods 
which are of high accuracy and efficiency because of their excellent approximation 
properties for smooth functions and the existence of fast algorithms for the Fourier 
transformation. The power law tails of the solution, as encountered in the preceding 
section, however will lead to an increasing number of "echoes" as time goes on, due 
to the chosen periodic domain of computation. For sufficiently large periods though 
these echoes will be small on the studied time-scales and they will not influence the 
analysis of the KP oscillations. 

We use here an adapted version of Trefethen's code for the KdV equation (Chap. 10 
in ^Uj) which is available at 2D- The basic idea of the code is the use of a discrete 
Fourier transform in the spatial coordinates as well as of an integrating factor 
such that the time derivative is the only linear term appearing in the equation. 
This preprocessing, as introduced in [18) . leads to an equation without a stiff part 
and allows for larger time steps. The method is slightly more efficient than the 
semi- implicit approach of J2|> or time splitting techniques such as 1TJ e.g.. More 
importantly, however, it allows for the use of higher order time discretizations which 
are convenient in the context of strong gradients, as needed in Subsection l5.1l below. 
Keep in mind that an n-th order approach introduces a numerical dispersion of the 
order (Ai) n+1 , which leads to a "pollution" of the high spatial frequencies in the 
numerical time evolution. On the other hand though we need these high frequencies 
to resolve numerically the solution near a gradient catastrophe in the solution of 
the dKP model. Thus it is helpful to be able to apply a higher order method in 
the time discretization within an efficient approach which allows to deal with 2 + 1 
dimensional settings. 

As before, let u(t) be the Fourier transform of u(t) and denote by u 2 the Fourier 
transform of u 2 . Then the Fourier transformed KP equation (|1.1|) reads 

(3.1) ^ + (^-^) 2 +5^ = °- 

Equation (|3.1|l has to be regularized for k x = in order to give numerically sense 
to the term —i/k x . Obviously division by zero is not allowed, even if the result 
would make sense analytically in a certain limiting procedure. To carry out such a 
limit numerically one has to use appropriately chosen finite numbers. Here this is 
done by adding to k x in the denominator a small imaginary part of appropriate sign 
(corresponding to A). In the numerics we add the smallest floating point number 
which Matlab can represent, 2.2 . . . 1CP 16 . Equation 13. 1|) is then equivalent to 

(3.2) d t ( e it ( Xk « /(k * +iX0) - t2k * ) d > J + — e it ( Afe ?/( fc -+ iA °)- e2fc x) ^2 = 

To solve equation Il-S.21 numerically we use the Fast Fourier Transform (FFT) in 
MATLAB for the dependence on the spatial coordinates and a fourth-order Runge- 
Kutta method for the time integration for the reasons given above. The important 
term in this integration is e - lAt (M y /{k x +i\o)-e k J as can be seen, e.g., from the 
simple time discretization 

(3.3) U(t + At) = c -iAt(A^/(fe,+iA0)-^) ^ {f) _ - 
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Since we use an explicit method for time integration, stability is an issue. We find 
that due to the integrating factor, time steps of the order At ~ l/(N x N y ) lead to 
a stable time evolution. In the following, the empirically found values for the time 
steps needed in the respective computation are always given. 

As already discussed in Section [5J we require our initial data to be subject to the 
constraint (|2.1U|) . In the periodic setting we analogously impose 

Ul (x,y)dx = 0, 

-jri B 

where 2itL x is the period in x. For such data we solve numerically for the function 
v defined by u = ik x v to reduce numerical errors. 

Remark 3.1. The numerical code is also able to propagate initial data which do 
not satisfy the constraint i|3.4|) . which however yields a non-smooth solution in 
arbitrary short times. The reason for this is the analytical behavior of the term 
e -itAfc H /(fe^+iAO) a pp earm g m (|3,3|) , Indeed, for k x — 0, this term is numerically 
equal to 0, unless k y = where it is equal to 1. This leads to a continuous but not 
diffcrcntiable solution at x = 0, see Fig.^ 



ia-. 

IB- 

-ja -. 

12- 

a . 




Figure 1. Solution at tim e t = 4.8 x 10~ 4 to the KP-I solution 
with initial data 1/ cosh 2 (-\/ x 2 + y 2 ) which are not subject to the 
constraint 

3.2. Numerical test cases on solitons. To test the accuracy of the code we 
compare, for the different time steps, the numerical solution to an explicitly known 
exact solution of the KP equation. Since the KP equation is completely integrable, 
large classes of explicit solutions are known. The most discussed cases certainly are 
the above-mentioned solitons. We use here the 1-soliton The propagation 

of this solution formally does not test the y-dependent terms in the KP equation. 
However, this is only partially true since the line solitons are known to be unstable 
against perturbations in the case of the KP-I equation. Since unavoidable numerical 
errors provide some form of perturbation, the propagation of l|2.6(l also tests the 
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y-dcpendent terms in the numerical implementation, see the example below. The 
reason for the use of a y-independent solution is that line solitons with an explicit 
y-dependence, cf. will not be periodic in y. 

The test of the propagation of the 1-soliton initial data is performed for the KP- 
I model in the following setting: The x-coordinate takes values in the interval 
[—irL x ,irL x ], where the length L x is always chosen such that the coefficients for 
the initial data are of the same order as the rounding error for high frequencies 1 . 
This reduces the error due to the non-periodicity of the initial data at the interval 
junctions to the order of the rounding error. The choice of the corresponding length 
L y is not important in this context, thus we shall choose it to be of the same value 
as L x . We use the 1-soliton solution at time t = with L x = 10 and xq = —L x in 
the initial data and determine for each time step the difference between the exact 
and the numerical solution. The computation is carried out with N x — 2 11 and 
N y = 2 7 modes (we will always use powers of 2 here since the FFT algorithm is 
most efficient in this case, but this is not necessary) and 4 x 10 3 time steps for 
t G [0,6]. The L°° norm of the difference between the numerical and the exact 
solution is shown as a function of time in Fig. 



-6 




1 2 3 4 5 6 

x 



Figure 2 . Numerical errors for the time evolution of 1-soliton ini- 
tial data: L°° norm of the difference between exact and numerical 
solution (maxdiff) and deviation from mass conservation (err). 

An alternative test of the numerical precision is provided by conserved quantities. 
Due to unavoidable numerical errors such quantities will be numerically a function 
of time. Typically the energy (|2.3H is a convenient quantity in this context. However 
the term d~ l d y u in general will have the same analytical properties as KP solutions 
not subject to the constraint Ij3.4|l . namely a cusp along the ir-axis. Since we 
use spectral methods, this non-analyticity of the integrand in (|2.3(l would lead to 
numerical errors, which however are mainly due to the way the integral in l|2.3|) is 
evaluated, and not to the numerical error in u(t) itself. 

^MATLAB works internally with a precision of 1CU 16 . Due to rounding errors machine preci- 
sion is generally limited to the order of 10 — 14 . 
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Therefore we shall instead consider the mass H2.2[) . which is also conserved in time, 
but which admits an integrand with better regularity properties than the energy. 
Since mass conservation is not implemented, it provides a strong test for the accu- 
racy of the code. We define the corresponding numerical error function by 

(3.5) -(fl-l-^f, 

where M[w(i)] is the numerically calculated mass which is obtained via FFT (in both 
spatial coordinates) of the integrand in (|2.2[1 at each time step. For the example 
of the 1-soliton, this function is shown in Fig. |5J It can be seen that the error 
obtained via the integral quantity is typically one order of magnitude higher than 
the maximal local difference of the exact and the numerical solution. In cases where 
no exact solution is known, we shall use mass conservation as an indicator of the 
precision of the numerics. Consequently, the number of modes and the time steps 
will always be chosen in a way such that the value of the function err(i) is at least 
one order of magnitude lower than the precision of the numerical solution we are 
aiming at. 

As an example for the instability of the above given line-soliton in the case of the 
KP-I equation, we consider, as in [33], a strongly perturbed line-soliton of 
with initial data 

12 

(3.6) u\(x,y) = J • 

cosh (x — xq + <5cos(0.2y)) 

In Fig. we show the growth of the perturbation with time and the formation of 
lump solitons for 5 = 0.4, L x = 12, L y = 10, At = 1.33 x 10~ 4 and N x = N y = 512. 
Due to the algebraic decay of the lump solitons for |a;|, \y\ — > oo, relative mass 
conservation is of the order of 10 -3,5 in this case. This is the precision with which 
the code propagates lump solitons for the given parameters. 

Remark 3.2. Notice that the formation of power-law tails in the time evolution of 
initial data with compact support will lead to a Gibbs phenomenon at the bound- 
aries in our spatially periodic setting. These effects are small for sufficiently large 
periods, but clearly noticeable. Thus in the examples studied here, relative mass 
conservation will be of the order of 10 -7 , even for very small times, whereas in 
the corresponding computations for the KdV equation it is of the order of machine 
precision [lfij . On the other hand, on the longer time scales considered here (when 
studying the rapid modulated oscillations we are interested in), it is still possible 
to achieve relative mass conservations of the order of 10 . The same holds for 
the propagation of lump solitons such as (|2.5fl . where the slow (algebraic) decay 
as \x\, \y\ — > oo again leads to Gibbs phenomena at the boundaries. Nonetheless 
our scheme is able to propagate lump solitons on standard computers with relative 
mass conservations of the order of 10~ 4 as the code in |12| . 



4. Modulation theory for small amplitude solutions 

Since the small dispersion limit e — > of i|l.lf> is only poorly understood, we first 
consider the analogous problem for solutions which are asymptotically small, i.e. of 
order 0(e). To this end we perform a (formal) multiple scales expansion, similar to 
those given in |3fil I4fi| . for such kind of solutions. Since the quadratic nonlinearity 
in this case is of lower order, we expect that the appearing fast oscillations can be 
described by a rather simple ansatz. More precisely they will be purely x-dependent. 
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Figure 3. Time evolution of a perturbed line-soliton, growth of 
the perturbation and formation of lump solitons in the KP-I equa- 
tion. 



Namely, via the so-called Airy equation with plane wave initial data, i.e. 

(4.1) dtv + d xxx v = 0, u| t =o = e ' We - 

Note that this yields an asymptotic description which does not involve fast oscilla- 
tions in the y- variable, cf. below. This is certainly a simplification which neverthe- 
less is justified by the computations given in Section [21 for the linear equation and 
by our numerical results. 

4.1. Multiple scales expansion. To be more precise, let us consider the scaled 
KP equation Ijl.ljl . assuming that its solutions admit an asymptotic expansion of 
the form 

(4.2) u{t,x,y) ~ eu (t,x,y) +e 2 Ul (t,x iy ) + 0(e 3 ), e < 1. 
As usual in multiple scales expansion methods, we consider 

(4.3) Uj = u j {j i t,x,y ) X ) Y,T) 1 V j e N, 

introducing again the fast variables X — x/e, Y — y /e, T — t/e, as well as the 
additional time-scale r = et. Moreover we assume Uj to be periodic w.r.t to X, Y, T. 
Note however that in the initial data (|4.1(l we do not take into account the fast scale 
Y = y/e. Thus we shall neglect this dependence in the following and one easily 
checks that this is consistent with our asymptotic expansion up to any order in e. 
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Also note that we have to take into account additionally the (slowly varying) time 
scale t = et in order to balance the nonlinearity with the term d x dtu. 

Equating powers of e in our expansion, we find that the leading order term uq is 
given by 

(4.4) u (t,x,y)=Met,x + 3r ] 2 t,y)e i ^ x+u ^ t y e + c.c., ^(77) = t? 3 , 

for any given 77 G R, 77 ^ 0, where ?/>o(i, x), y) G C and "c.c." refers to the 
complex conjugate. The fast oscillations are henceforth described by plane waves 
propagating according to io{rf) = rj , the dispersion relation of the Airy equation. 
Moreover, in (|4.4ll . we need to take into account a drift defined by 

£(t, x) := x + 3?7 2 i, 

i.e. we need to perform a change of variables into a reference frame moving at 
group velocity v g — 3?y 2 . Thus, on the so-called ballistic scales (t,x) = (eT, eX), 
the slowly varying amplitude satisfies the transport equation 

(4.5) d t TPo ~ 3t7 2 c^Vo = 0, 

which is obtained from the solvability condition for terms of order 0(e). 

Consequently one finds that the fast scale dynamics in the first order corrector u\ 
is determined via 

(4.6) d x (d T + d xxx ) Ul = -\d xx {i> 2 (r^(t,x),y)e 2i ^ x+u ^ +cc.) 
and thus u\ is given by 

Ul (t,x,y) = -L^ ( Tj Z(t,x),y)e 2i ^+^)T) + M T^(t,x),y)e^ x+ ^)T) 

(4.7) 67] 

+ c.c. + (j>(T,£(t,x),y) 

T =et,X=x/e,T=t/e 

where <f> can now be interpreted as a mean field generated by ipQ. From the mathe- 
matical point of view it is necessary to take into account this non-oscillatory field to 
balance the terms generated via the quadratic nonlinearity. Wc eventually find that 
ipQ and 4> solve a coupled system of Davey-Stewartson type (DS). More precisely we 
get 

J id T ip ~ 3??%V>o + - dyyipo - ( \ip \ 2 + r\4> ) "00 = 0, 

(4.8) < V V677 J 

[ 3»7 2 %0 + A9 w </. + %|Vo| 2 = O. 

This mean field system describes the dispersive effects which are visible on time- 
scales of order 0(et). Analogous to the KP equation, the DS system represents 
a completely integrable model in d = 2 spatial dimensions. The case where, for 
i] > 0, A = +1, is referred to as the "hyperbolic-elliptic" case, cf. whereas 
for A = —1 we are in the so-called "elliptic- hyperbolic" situation. Here we shall 
only focus on the former case, induced by the KP-II model, as the latter requires a 
special numerical treatment due to the appearing wave-type operator in the second 
equation of (|4.8(l . Thus we require, that the mean field behaves like 

(4.9) #T,£,y)->0, as \y\ -> 00. 

Remark 4.1. An analogous asymptotic expansion for the KdV equation yields 
a cubic nonlinear Schrodingcr equation, instead of the DS system, as has been 
rigorously proved in |35) . We expect that a similar analytical strategy could also 
be applied in our case to rigorously establish the above given formal asymptotics. 
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In the present work, the numerical algorithm used to solve 1)4. 8[l is analogous to the 
one used for the KP equation (for an alternative approach see @]). More precisely, 

denote by i/jq(t, %j ky) the Fourier transform of ipo(r, £, y) and let V'olV'ol an d "00^ 
be the Fourier transforms of V'olV'ol anc ^ ^Po^i respectively. Thus we get for the 
system 1)4. 8|l 



(4.10) 



ml 



with 77 > 0. Again we shall use an integrating factor to avoid a stiff part in the first 
equation of l|4.10|l . i.e. 

(4.11) d T (e^ 3 "^ 2 -^/'"^) + ie- iT ( 3 " fc2 - Afe » Q-^f| + rj^j = 0. 

The time integration will be carried out as before with a fourth order Runge-Kutta 
method. Accuracy of the code will be checked by conservation of the wave energy 

(4.12) N^o(r)] := / \Mr,^y)\ 2 d£dy, 
which is preserved in time for the DS system 14.8fl . 

4.2. Numerical examples. To study a concrete example we consider, for rj = 1, 
real-valued initial data to the DS system (|4.8I) in the following form 



(4.13) iJj \ t=0 = M^,y)^-dxSech 2 (R), R := ^x^+y 2 . 

The choice for these initial data is motivated by the fact that they are smooth, 
localized in x, y, and that they could, in principle, also serve as initial data for the 
KP equation 1)1.1)1 . since they satisfy the constraint (|3.4|) . 

We then solve the DS system (|4.8ll with initial data (|4.13|) for t e [0, 0.4] with N x = 
N y = 512 and At = 2 x 10~ 3 . The wave energy is conserved in this computation 
up to the order of machine precision (the error is smaller than 10 -13 ). We show 
the real part of the function ijjo(t) for several values of t in Fig.0] It can be seen 
that the initially localized pulse spreads and exhibits oscillations, both mainly in 
x-direction. For longer times, the modulation of the DS solution becomes more 
pronounced as can be seen in Fig. [SJ Notice that the real and imaginary part of 
ipo(t) have almost the same envelope, they mainly differ by a phase shift. This 
is more obvious by considering the absolute value of ipg in Fig. which virtually 
shows no modulations. 

Corresponding to (|4.13l) , the KP-II equation is then solved with initial data 

(4.14) Ul (x,y) = 2e ^{Rjcos (^) , R=^x 2 +y 2 . 

Note however, that these initial data do not satisfy the constraint (|3.4|l . To enforce 
(13.41) we numerically compute the Fourier transform of u\ and set the Fourier coeffi- 
cients of all terms corresponding to k x = in these data equal to zero. This can be 
justified by the fact that in our examples, where we take e = 0.1 and smaller, these 
Fourier coefficients are only of the order 10 -14 , i.e. of the order of the rounding er- 
ror. Thus, within numerical precision, we can directly use these initial data, which 
are now adapted to both, the asymptotic expansion and the integral constraint 
(13. 4|) . For the precise parameters used in the computations of the KP equation, for 
times between and 4, see table 
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Figure 4. Time evolution of Reipo(t), solution to the DS system 
with initial data l|4.13|) . shown for several values of t > 0. 
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Figure 5. Keipo(t) an d Im V'o(i), obtained from the initial data 
(|4~T3|) . at time t = 2. 
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Figure 6. | Vo (*) 1 2 , obtained from the initial data Q4.13|) . at time 
t = 2. 
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Table 1. Parameters in the computation of the KP solutions 



In Fig. U| we plot the solution of the KP equation for e = 0.1 for several values 
of t. As expected the pulse travels roughly with the velocity v g — 3 (notice that 
the shown region in the plots is co-moving as can be seen from the x-scales) and 
the initially localized pulse spreads mainly in x-direction. Since we only consider 
small amplitudes here, the characteristic KP tails are hardly visible on the studied 
timescales. Thus the periodic boundary conditions do not influence the model. 

The above initial data are of the required form for the asymptotic expansion (|4.2|) 
in leading order of e. In principle one could try to obtain a refined asymptotics, 
including higher order corrector terms. This however yields the highly nontrivial 
problem of satisfying the KP constraint l|3.4|l up to sufficient high order. More 
precisely if one takes into account also the first order corrector Ij4.7|l . then, since 
4>{t) is proportional to the Fourier transform of I'i/'oWl 2 , the coefficients for k x — 
are of the order 0(e 2 ). In order to satisfy the constraint also up to 0(e 2 ) errors, 
it consequently would be necessary either to choose ip2(0,x,y) in an appropriate 
way, or to consider a function tpx(x,y) which is subject to a nonlinear integral 
constraint following from 1)3. 4[l . In the following we shall simply neglect higher 
order corrector terms in the initial data as we only aim (numerically) for the leading 
order asymptotic description. 
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Figure 7. Time evolution of the initial data (|4.14|> . governed by 
the KP-II equation with e = 0.1, shown for several values of t > 
and shifted by x — * x — 12. 



The above given expansion indicates that the solution to the DS system l|4.8[L 
with initial data ipi, should approximate, via (|4.2|) and (|4.4|> . the corresponding KP 
solution up to errors of order 0(e a ), for some a > 1, on any finite time-interval. In 
Fig. [HI we plot the approximating solution 

(4.15) u app {t,x,y) = 2eRe (V(ei, £(t, x), V )e l{x+t)/c ) , 

for several values of t. It can be seen that this asymptotic solution gives the expected 
description of the KP solution, which is even better visible in Fig. El where we have 
plotted the two solutions for t = 4 and y — in one frame. The difference of these 
two solutions for t = 4 can be seen for y — in Fig. 1101 and in the whole (x, y)-plane 
in Fig. El 

Remark 4.2. Notice the difference between the KP solution and the asymptotic 
solution for different signs of the co-moving coordinate £(t, x). We chose initial data 
which are odd in x. Obviously the KP equation does not conserve the parity of the 
initial data, whereas the DS system does. In our example the asymptotic solution 
is, up to the leading order in e, an odd function in x, whereas this is no longer 
true for the KP solution. Higher order terms in the asymptotic solution, which we 
neglected here, will break this symmetry. 

For smaller values of e the approximation gets better, but the numerical resolution 
of the high frequencies becomes increasingly more difficult. The main problem is 
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Figure 8. The asymptotic solution u app (t,x,y) = 
2e Re (i/j(et, £(t, x), y)e 1 ( x+ *)' e ) approximating the true KP-II 
solution of Fig. [7| shown for several values of t > 0. 

that the DS system approximates the KP solution only on the time scale et. Thus 
in order to see the effects of the modulations due to the DS system within the KP 
solution one would need to solve the KP equation for extremely long times if e is 
chosen very small. This is of course numerically rather expensive, in particular in 
our 2 + 1 dimensional case. We will therefore limit our analysis to the range 0.01 < 
e < 0.1. The above given example however clearly illustrates the applicability of 
the asymptotic expansion. For e = 0.01 the point wise difference between the KP 
solution and the asymptotic solution is, as expected, at most of the order 10 and 
thus barely one order of magnitude higher than the numerical error, cf. Fig. 1121 To 
get more insight we plot the L°°(M. 2 ) norm of the difference, denoted by A 00 (t), in 
dependence of e in Fig. 1131 This error is found to be monotonically increasing in 
time and furthermore it decreases roughly like e 9 / 4 . Indeed the data obtained at the 
last time-step t = 4, can be fitted by a straight line — log 10 Aqc = —a log 10 e + b with 
a = 2.27 and b = —0.58. The correlation coefficient is then found to be r = 0.999, 
the standard error for a is a a — 0.08. 

To define an integral quantity A2 (t) for the error between the true KP solution and 
its (leading order) asymptotic description (|4.15(l . we also consider the L 2 (R 2 ) norm 
of the difference. Indeed this is the most widely used definition of an asymptotic 
error in such singular limiting regimes, cf. |35| for the analogous definition in the 
KdV case. Thus we integrate the square of the difference of these solutions over the 
domain of definition via an FFT (in both variables) and normalize this quantity 
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FIGURE 9. Solution to the KP-II equation with initial data (|4.14|) 
(blue) and the corresponding asymptotic solution (green) for y = 0, 
t = 4 , and e = 0.1. 
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Figure 10. Difference of the solution to the KP-II equation with 
initial data (|4.14() and the corresponding asymptotic solution for 
y = 0, t = 4, and e = 0.1. 

dividing it by the area of the fundamental domain, i.e. 

(4.16) A 2 {t) := 7== ( / / (u(t,x,y) - u am (t,x,y)f &x &y 

Z-KyJ L x L y \J-nL x J-nL y 
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Figure 1 1 . Difference of the solution to the KP-II equation with 
initial data (|4.14|) and the corresponding asymptotic solution for 
t = 4 and e = 0.1. 




Figure 12. Difference of the solution to the KP-II equation with 
initial data (|4.14(l and the corresponding asymptotic solution for 
t = 4, and e = 1(T 2 . 



We find that A2(t) again increases monotonically in time. For t = 4 we get the 
values of A2 shown in Fig. El By linear regression the data can again be fitted 
by a straight line — log 10 A2 = — alog 10 e + b, with a = 2.48 and b — 0.53. The 
correlation coefficient in this case is r = 0.997 and the standard error for a is 
a a = 0.17. 
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Figure 13. Error Aoo for several values of e. The data can be 
fitted by least square analysis with a straight line — log 10 Aoo = 
— alog 10 e + b with a = 2.27 and b = —0.58. The correlation 
coefficient is r = 0.999, the standard error for a is a a = 0.08. 

Remark 4.3. It is not surprising that A2 decreases faster with e than A m since 
not only do the amplitudes shrink with e, but also the pulse is more localized in x 
and y. Thus a smaller integral error is to be expected. 

In other words we get that numerically the L 2 (M 2 ) difference between the true KP 
solution and its asymptotic description is approximately of the order 0(e 5 / 2 ). This 
is remarkable since it fits with the analytical results of where the analogous 
limit from KdV to the cubic nonlinear Schrodinger equation is considered. 

Remark 4.4. Note however, that for a comparison of our asymptotic description 
(|4.4|l with the one given in [23 one has to take into account a rescaling of the 
spatial scales on which the modulation amplitudes varies. Thus rescaling x and y 
(the latter scale is of course not present in one has to plug this in "by hand") 
such that the two descriptions match each other, one realizes that an additional 
factor e has to be taken into account in order to conserve the L 2 (K 2 ) norm of the 
solutions. In summary one checks that Theorem 1 in |ri5| yields an asymptotic error 
of the order 0(e 5 / 2 ), i.e. exactly as in our case. 

5. Small dispersion KP solutions with amplitudes of order 0(1) 

In this section we will study numerically solutions to the KP equation in the regime 
of small dispersion e <C 1. In contrast to Section 01 though we shall not restrict 
ourselves to small amplitude solutions. Rather we shall consider smooth initial data 
with an amplitude of order 0(1), as e — > 0. To this end let us first discuss briefly 
the resulting (formal) limiting equation (|1.8f) and its numerical implementation in 
the next subsection. We shall then investigate a concrete example with initial data 
of the form 

(5.1) uifotf) = -6d x sech 2 (R), R = ^x 2 +y 2 , 
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Figure 14. Error A2 as defined in 14.161 for several values of e. 
The data can be fitted by least square analysis with a straight 
line — log 10 A2 = — alog 10 e + b with a = 2.48 and b = 0.53. The 
correlation coefficient is r = 0.997, the standard error for a is 
a a = 0.17. 

as can be seen in Fig. El Note that these data satisfy the constraint p.lOJl ■ The 
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Figure 15. Initial conditions for the KP equation. 



factor 6 in H5.1(l is included in order to be able to compare our results with those 
given in ^1 f° r the KdV equation, where a different definition of the amplitude 
u(t) was used (the u(t) there corresponds to 6u(t) here). With the introduction of 
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the factor 6 we expect the same dynamical time scales in the KdV sector of 
as in [H>| . 

5.1. The dKP equation and its dissipative regularization. We do not expect 

the formal limiting dKP equation (jl.8|l . i.e. 

d x (dtu + u d x u) + A dyyU = 0, 

to be the correct description of the dispersionless limit for (O) , at least not for all 
times t > 0. This believe stems from the closely related situation encountered in 
the inviscid Burgers, or Hopf equation, given by 

(5.2) dtu + u d x u = 0, u\ t _ = u\{x). 

Equation (|5.2[1 can be seen as the (formal) dispersionless limit of the KdV equation. 
In general though this only holds for some finite time < t c < oo. More precisely, 
the description of the disperionless KdV limit via (|5.2|1 fails, after the appearance 
of the first shock in the solution of (|5.2(l . The corresponding break time is given by 



(5.3) t c = min 



io£S\ d x ui(xo), 

which can be easily seen when solving l|5.2|l by the method of characteristics, yield- 
ing 

(5.4) u(t, x) — ui(x ), x — ui(x )t + x q , 

the so-called Hopf solution. Consequently, we also expect in the case of the dKP 
equation that a shock will be formed after some finite time. (Clearly the behavior 
of the dKP equation is completely analogous to the Burger's case, if one considers 
y-independent solutions.) On the other hand, for small, but still finite, e <^ 1, this 
"dKP-shock" presumably will be smoothed out by rapid oscillations, similar to the 
KdV case, cf. the numerical examples below. 

Remark 5.1. To convince ourselves that, apart from the KdV sector, the concept 
of shocks is really applicable in the dKP equation we remark that its principal part 
is given by 

(5.5) Pu := d xt U -\- 11 XX U -\- XOyyU. 

As one easily checks, the coefficient matrix in this partial differential operator P 
admits three, distinct, real eigenvalues given by 

(5.6) /xi )2 = ^(u± \Jv? + l) , ^ 3 = A. 

Thus the dKP equation indeed falls into the class of second order hyperbolic PDEs.c 
However the associated initial value problem as considered in this work, furnishes 
a characteristic Cauchy problem, a fact which has already been discussed in Sub- 
section 

To circumvent the problem of shock solutions in the formal limiting dKP model 
we shall use a dissipative regularization. More precisely we add a small dissipative 
term, in the form 

(5.7) d x (dtu + u d x u — a d xx u) + A d yy u = 0, A = ±1. 

with < a <g; 1 being some small real- valued parameter, such that a ~ 0(1), as 
e — > 0. Equation (|5.7|l can now be seen as the KP analog of the viscous Burgers 
equation, i.e. 

(5.8) dtu + u d x u — a d xx u = 0. 
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As in %£2fy we obtain that l|5.7|l can be solved via 

(5.9) d t (eW/^+'^+^u) + 1 -k x e *P^/(*-+ i ^)+^) ^ = 0. 

Obviously equation l|5.7|l is no longer conservative, but the dissipative term will 
smooth out the shocks of the dKP equation. In other words, as a — > 0, we expect 
the solution of (|5.7|) to tend to some kind of entropy solution. 

Remark 5.2. Without the dissipative regularization, the numerical code would 
break down shortly before the formation of shocks, due to the appearing steep 
gradients which cause instabilities in our explicit time integration scheme. However, 
even with an implicit, unconditionally stable method, such as Crank-Nicholson, the 
numerical code would break down close to a shock. The reason for this is mainly 
due to the so-called aliasing error, i.e. a pollution of the spectral coefficients by 
high frequencies, see, e.g., HJ. Because of the nonlinearity in the KP equation, this 
high frequency noise yields severe problems near the gradient catastrophe. Even a 
de-aliasing, as used in, e.g., |16| for the KdV case, will not stabilize the code close 
to the breakup, since the aliasing error simply cannot be suppressed there. 

Since the dKP equation is also completely integrable, many explicit solutions are 
known, cf. ^3 [213 [57]. m our study, though, we are interested in rather general 
initial data, where no closed form of the corresponding solutions is known. In 
the case ()5.1JI the corresponding KP-I solution (A = —1) is shown in Fig. 1161 
where a = 0.01. The computation was carried out with N x = 4096, N y = 128, 




Figure 16. Solution to the regularized dKP-I equation with ini- 
tial data l|5.1(l and a = 0.01, plotted at time t = 0.3. 

L x = L y = 10 and At — 5 x 10~ 5 . Due to the dissipative term cx ad xx u the mass 
is no longer conserved, but in the given example the loss is only of the order of a 
few percent. Thus the shown solution should indeed be very close to a true shock 
solution of the dKP equation. Note that the tails, as x — > oo, are clearly visible. 
Out focus here is now on the the region of steep gradients in the solutions to dKP. 
To this end it can be seen from Fig. 1171 that the derivative d x u is always maximal 
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Figure 17. The x-derivative of the regularized dKP-I solution 
with initial data l|5.1|l and a = 0.01, plotted at the time t = 0.3. 

on the a:-axis (note that we plot here —d x u). As expected the maximum of the 
derivative is of the order 1/cx, and it is significantly bigger in the wave front for 
x > 0, where the tails form. This can also be inferred from Fig. 1181 where the 
dKP-I solution is plotted on the x-axis for several values of the time. We observe 

t= 0.000 t= 0.225 




-4 -2 2 4 -4 -2 2 4 




-4> • • • 1 -4 1 • • • 1 

-4 -2 2 4 -4 -2 2 4 



Figure 18. Solution to the regularized dKP-I equation with ini- 
tial data (|5.1() and a — 0.01, plotted on the x-axis for several values 
oit. 

that the gradient catastrophe is reached for x > roughly at the time t c « 0.23. 
For t > t c , the gradient remains more or less unchanged due to the implemented 
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dissipation in (|5.7(l . The gradient catastrophe for the second wave front for x < 
is reached roughly for t ss 0.3. 

For A = 1, i.e. the KP-II case, the situation is similar, but the tails are now 
directed towards x — > — oo. The gradient catastrophe is now reached first in the 
wave front for x < 0. This is due to the fact that the dKP equation ljl.8|l . as well 
as its dissipative regularization l|5.7|l . are invariant under the simultaneous change 
of x — > — x, u(t) — > — u(t), and A — > —A, for t > 0. For initial data of the type 
(|5.1|) which are odd in x, i.e. uj(—x,y) = —ui(x,y), one consequently obtains the 
dKP-II solution u + (t, x, y) from the dKP-I solution u_(t, x, y), with the same initial 
data, by identifying u + {t,x,y) = —U-(t,—x,y). 

5.2. KP oscillations. In the following we shall study the oscillations for small, 
but still non-zero, e <C 1 in the solution of ll.lfl . 

In the closely related KdV case this limiting regime is rather well understood, see, 
e.g., ^JEHlEniEII] : Roughly speaking, the corresponding analytical approach con- 
sists of the following steps: After the breaking in the corresponding Hopf solution 
(|5.4() . an oscillatory zone is identified via the solution of the so-called Whitham 
equations, which depend only on the slow coordinates x and t. The corresponding 
approximate solution to the KdV equation in this region is consequently given in 
terms of theta functions, evaluated on the fast coordinates X — x/e, T = t/e, 
whereas the branch points of the underlying Riemann surface are determined via 
Whitham's modulation equations. Outside the oscillatory zone, the KdV solution 
is approximated by the corresponding solution of the Hopf equation. For a com- 
parison of this asymptotic solution with a numerical KdV solution see |lfi| . 

There are several obstacles to overcome when one tries an analogous approach in the 
KP case. First, one should keep in mind that the formal dispersionless KdV model, 
i.e. the Burger's equation (|5.2|) . is a first order PDE which can be integrated by the 
method of characteristics. On the other hand, the corresponding dKP equation, is 
a second order PDE. Although the dKP equation is completely integrable and thus 
explicit solutions are known, there exists no general procedure of integrating this 
equation for generic initial data so far. Secondly, the Whitham equations in the 
KdV case can be brought, via the hodograph transform, into the so-called Riemann 
linear form. They can then be again solved (at least implicitly) for generic initial 
data. Thus all quantities entering the theta functional solution, as for instance 
the phase, are known explicitly via their dependence on the branch points of the 
Riemann surface defined by the solution of the Whitham equations. 

Remark 5.3. In the KP case, the corresponding Whitham equations were solved 
formally in [28] . So far, however, it is not clear how to bring this solution into a 
form which is convenient for a numerical simulation. This is especially true for the 
non-monotonous initial data with two inflection points we are considering here, a 
case which has not even been studied in the much better understood KdV setting. 
(Note that the two inflection points are needed to satisfy the constraint Q2.10JI .) In 
addition it would be necessary to have an explicit formula for the corresponding 
phase entering the theta function to be able to compare the asymptotic solution to 
the exact solution in the Whitham zone. 

We hope that the numerical results presented in this paper help to develop a similar 
asymptotic description as in the KdV case which will be the subject of further 
research. 
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Now, we shall first show (numerically) that the solution of the dKP equation gives 
the correct limiting behavior as e — > before breakup. To this end we compare 
the dKP-I solution UdKp(i), with initial data (|5.1|l . and the corresponding KP-I 
solution UKp(t) at times t < 0.2, i.e. up to times close to (but still before) the 
appearance of the first shock in the dKP-I solution. We shall consider values of e 
between 0.1 and 0.01, with N x = 2 U , N y = 2 7 , and At = 2 x 10 -5 . Note that we 
solve here the dKP equation without dissipative regularization since we stop the 
computation before the appearance of the gradient catastrophe. We also remark 
that this is possible with a relative mass conservation of 10~ 4 39 . 

We find that the difference between the solutions UdKp(^), t < t c , and ukp(0 
increases monotonically with time. In the considered range of e and for t = 0.2, 
the L 2 norm difference A2, defined in (|4.16|l . decreases roughly like 0(e 3 / 2 ), as 
e — » 0. More precisely we find that the data can be fitted by a straight line 
— log 10 A2 = — alog 10 e + 6, with a = 1.45 and b — 1.30, as can be seen in Fig. ED 
The correlation coefficient is found to be r = 0.999, the standard error for a is 
o a = 0.056. 

4.5 1 1 1 1 1 1 1 1 1 1 1 




2.5 1 1 1 1 1 1 1 1 1 1 1 

1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 

-lofl 10 e 

Figure 19. Error A2 for several values of e. The data can be 
fitted by least square analysis with a straight line — log 10 A2 = 
— alog 10 e + b with a = 1.45 and b = 1.30. The correlation coeffi- 
cient is r = 0.999, the standard error for a is a a = 0.056. 

Remark 5.4. Note that the 0(e 3 / 2 ) error in A2 fits with our results of Section^] 
There we numerically found that the A2 error for asymptotically small solutions, 
i.e. solutions of the order 0(e), decays like 0(e 5//2 ) for small Kl. 

For the corresponding L°° norm difference Aoo we find numerically that the error 
roughly decreases like e. More precisely again, we get (see Fig. |2"0"|) . that the data 
can be fitted by a straight line — log 10 A2 = — alog 10 e + b with a = .96 and 
b = —0.31. The correlation coefficient is found to be r = 0.993, the standard error 
for a is a a = 0.099. Notice that at t = 0.2, the first oscillations of the KP solutions 
appear as can be seen below. This leads to considerably large values of Aoo since 
the solution to the dKP equation will never show oscillations. 
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Figure 20. Error Aoo for several values of e. The data can be 
fitted by least square analysis with a straight line — log 10 Aoo = 
— alog 10 e + b with a = 0.96 and b = — 0.31. The correlation 
coefficient is r — 0.993, the standard error for a is a a = 0.099. 

For times t > t c , i.e. beyond the onset-time of the first gradient catastrophe, 
we shall study the KP-I solution, corresponding to the initial condition (|5.1f> . for 
0.01 < e < 0.1 and times t < 0.4. The computations are carried out with N x = 2048 
to N x = 8192, N y = 128, L x = L y = 5, and At = 2 x 10" 5 to At = 4 x 10~ 6 . 
Thereby we observe a relative mass conservation of the order of 10~ 4 to 10~ 3 . The 
following observations are now in order: 

• It can be seen in Fig. for the KP-I case that oscillations mainly form in 
the region of large (spatial) gradients appearing in the dKP solution given 
above. 

• As in the KdV case the number of oscillations increases with decreasing e, 
whereas the wavelength decreases, see Fig. l2"2l 

• As expected there are two oscillatory regions, one for x < and one for 
x > 0, both of which are shown in detail in Fig. 1231 To get a complementary 
view of the rapid oscillations in this region a contour plot of the same 
situation is shown in Fig. l2~4l 

• The oscillations are always at most rapid on the x-axis, see Fig. l25l and the 
results of the next subsection. 

We note that the fast oscillations are numerically well resolved, as can be inferred 
from Fig. |5HJ which shows two enlarged sections of the oscillatory zone. Moreover 
we observe that the time evolution of the initial data (|5.1|l is as expected from the 
analysis of the regularized dKP equation above. 

To obtain more insight, we shall focus in Fig. [53 on the KP-I solution plotted at 
y = 0, i.e. where the rapid oscillations appear first and where they are the most 
pronounced. We note that the first oscillation appears at the time t w 0.22, for 
positive values of x G K, i.e. shortly before reaching the gradient catastrophe of 
the dKP equation. On the other hand at time t rj 0.3, the first oscillation appears 
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Figure 2 1 . Solution of the KP-I equation obtained from the ini- 
tial data l|5.1|) for e = 0.1, at time t = 0.4. 
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Figure 22. Solution of the KP-I equation obtained from the ini- 
tial data l|5.1jl for e = 0.01, at time t = 0.4. 



at a negative value of the re-axis. Again this is shortly before a shock in the second 
wave front is reached. Both regions then develop more and more oscillations as 
time goes on until finally the situation described above and shown in Fig. 1251 and 
Fig. ESI is reached. In Fig.l2*31and Fig.|^onc can also see the y-dependence of the 
observed oscillations. One recognizes that although they are mainly concentrated 
on the 2-axis, there are some oscillations for small < \y\ <C 1 (see also the results 
of the next subsection) . 
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Figure 23. Solution of the KP-I equation obtained from the ini- 
tial data l|5.1|) for e = 0.01, at time t = 0.4. 



Remark 5.5. In the corresponding situation for the KdV equation, the oscillatory 
zone shrinks with e. The same is the case for the KP model as can be inferred from 
Fig. 1281 and Fig. [2SI Note that this shrinking is also present in the y-direction. 

In contrast to the dKP equation, the KP equation is not invariant under the trans- 
formation x — » —x, u(t) — > — u(t), and A — > —A. Thus the rapid oscillations for 
small e <C 1 in the KP-II equation will be rather different from the KP-I 
can be seen from Fig.OOl Here, the tails directed to x — > — oo are clearly visible as 
well as the strong oscillations for x < 0. On the other hand there are virtually no 
oscillations for x > 0. For the sake of comparison we plot in Fig. 1311 the situation 
for both KP models, i.e. for the cases A = ±1, at y = and t = 0.4. Whereas 
there are rather strong oscillations near both breaking wave fronts in the KP-I case 
(A = —1), the KP-II case (A = 1) enhances oscillations for x < while for for x > 
they are more or less suppressed. 

5.3. Dependence on the y-coordinate. As already mentioned, the KP equation 
was originally derived to describe quasi one-dimensional (dispersive) wave phenom- 
ena with only weak transverse effects. Thus physically interesting solutions will 
have a rather weak dependence on the coordinate y € R. Mathematically it is 
nevertheless interesting to study the effects of a stronger y-dependence. In this 
subsection we will thus have a closer look on this dependence in the context of 
low-dispersion. 
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Figure 24. Solution of the KP-I equation obtained from the ini- 
tial data l|5.1jl for e = 0.01, at time t = 0.4 
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Figure 25. Solution of the KP-I equation obtained from the ini- 
tial data l|5.1|) for e = 0.01, at time t = 0.4 and y = 0. 
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Figure 26. Solution to the KP-I equation obtained from the ini 
tial data l|5.1|) for e = 0.01, at time t = 0.4 and y = 0. 
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Figure 27. Solution of the KP-I equation obtained from the ini 
tial data H5.1(l for e = 0.01 for several values of the time. 
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Figure 28. Solution of the KP-I equation obtained from the ini- 
tial data l|5.1(l at time t = 0.4 for several values of e. 

Recalling the analysis of the Green's function for the linear part of the KP equation, 
given in Section 2.1, it appears that only the Airy-part leads to oscillations, whereas 
the term oc d yy u is responsible for the formation of tails. On the other hand, the 
gradient catastrophe in the dKP model is the result of the nonlinear term oc ud x u, 
as in the case of the Hopf equation. Consequently one expects for initial data with 
a non-trivial y-dependence that oscillations can only appear in a small vicinity of 
the x-axis, i.e. before the formation of tails takes over. This behavior has already 
been observed in the subsection above and it is further supported by Fig. 1321 which 
shows the KP solution for e = 0.01 and t = 0.4 in the vicinity of the x-axis. To 
study the effect of the dependence on the y-coordinate in more detail, we consider 
again initial data of the form 1)5.1(1 . However we shall now evaluate them at R = i?„, 
given by 

(5.10) Rl =x 2 + vy 2 , 

where now v £ is a deformation parameter. Obviously, for v = 0, there is no 
y-dependence, and the KP model reduces to the KdV equation. The parameter v 
thus allows for a continuous deformation of the KdV sector. 

In the following we only consider the case e = 0.1. The corresponding solution to the 
KdV equation can be seen in Fig. 123 For simplicity, we shall only compare solutions 
for different values of v plotted at y = 0, i.e. the region where the oscillations are 
the most pronounced. In Fig.|M|the solution to the KP-I equation with initial data 
(|5.1|) . (|5.10|) . is shown for several values of v. It can be seen that the case v = 0.01 




Figure 30. Solution of the KP-II equation obtained from the ini- 
tial data l|5.1|) for e = 0.1 at time t — 0.4. 
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Figure 31. Solutions to both KP models obtained from the same 
initial data H5.1(l . shown at time t — 0.4 and y = 0. 




Figure 32. Solution of the KP-I equation obtained from the ini- 
tial data l|5.1|) for e = 0.01 at time t = 0.4. 

is very close to the KdV situation. Thus a continuous deformation of the KdV 
sector to a nontrivial y-dependence is possible, also numerically. Increasing v has 
two effects: First we observe that the oscillations for x < are suppressed whereas 
the oscillations for x > are enhanced. This qualitative behavior has already 
been encountered above by comparing solutions of the KP-I and KP-II equation 
for the same initial data. Secondly, since a stronger y-dependence enhances the the 
formation of tails we observe that, the bigger v is, the more mass is transferred 
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Figure 33. Solution of the KdV equation with e = 0.1 obtained 
from the initial data (|5.1(1 with (|5.10|) and v — for the time 
t = 0.4. 

from the initial wave pulse to these tails. This implies that less mass is available 
to form shocks in the dispersionless case and thus rapid oscillations for small e are 
severely damped for large v, as is clearly visible in the case v = 10. 

Remark 5.6. In the case v = 10 the periodic boundary conditions produces very 
strong echoes, as it is clearly visible in Fig. |2U This however is the only example 
considered here, where the echoes dominate the oscillations which we want to study. 

Finally the corresponding transition in v for the KP-II case (A = 1) starting from 
a situation close to KdV is demonstrated in Fig. [35] There the oscillations for 
negative x are enhanced whereas the oscillations for positive x are suppressed as v 
increases. For large values of v we again find that more and more mass is transferred 
to the tails, i.e. pushed to x — > — oo. In summary we find that the tails tend to 
suppress strong gradients, and thus shocks in the dKP equation, as well as the 
corresponding rapid oscillations in (|l.lfl for small e. 

References 

1. M. J. Ablowitz and P. A. Clarkson, Solitons, nonlinear evolution equations and inverse 
scattering, Cambridge Univ. Press, Cambridge (1991). 

2. M. J. Abramowiz and I. A. Stegun (Eds.), Handbook of Mathematical Functions with For- 
mulas, Graphs, and Mathematical Tables, New York Dover 1972. 

3. J. C. Alexander, R. L. Pego, and R. L. Sachs, On the transverse instability of solitary waves 
in the Kadomtsev-Petviashvili equation, Phys. Lett. A 226 (1997), 187-192. 

4. C. Besse, N. Mauser, and H. P. Stimming, Numerical study of the Davey-Stewartson system, 
Math. Model. Numer. Anal. 38 (2004), no. 6, 1035-1054. 

5. V. N. Bogaevskii, On Korteweg-de Vries, Kadomtsev-Petviashvili, and Boussinesq equa- 
tions in the theory of modulations, Zh. Vychisl. Mat. i Mat. Fiz. 30 (1990), no. 10, 1487- 
1501; translation in U.S.S.R. Comput. Math, and Math. Phys. 30 (1991), no. 5, 148-159. 

6. M. Boiti, F. Pempinelli, and A. Pogrebkov, Properties of solutions of the Kadomtsev- 
Petviashvili I equation, J. Math. Phys. 35 (1994), issue 9, 4683-4718. 



3(5 



C. KLEIN, C. SPARBER, AND P. MARKOWICH 



v = 0.01 v = 0.1 




-6 -4 -2 2 4 -6 -4 -2 2 4 




-6 -4 -2 2 4 -6 -4 -2 2 4 



x x 



Figure 34. Solution of the KP-I equation obtained from the ini- 
tial data l|5.1(l . H5.1()|l for several values of v at time t = 0.4. 

7. A. dc Bernard and Y. Martel, Nonexistence of L 2 -compact solutions of the Kadomtsev- 
Petviashvili II equation, Math. Annalen. 328 (2004), 525-544. 

8. J. Bourgain, On the Cauchy problem for the Kadomtsev-Petviashvili equation, Geom. Funct. 
Anal. 3 (1993), no.4, 315-341. 

9. C. Canuto, M. Y. Hussaini, and A. Quarteroni, and T. A. Zang, Spectral Methods in Fluid 
Dynamics, Springer- Verlag, Berlin 1988. 

10. J. Colliander, C. Kenig, and G. Staffilani, Low regularity solutions for the Kadomtsev- 
Petviashvili I equation, Geom. Funct. Anal. 13 (2003), no.4, 737-794. 

11. T. Driscoll, A Composite RungeKutta Method for the Spectral Solution of Semilinear PDEs, 
J. Comput. Phys. 182 (2002), issue 2, 357-367. 

12. B. F. Feng, T. Kawahara, and T. A. Mitsui, conservative spectral method for several two- 
dimensional nonlinear wave equations, J. Comput. Phys. 153 (1999), no. 2, 467—487. 

13. E. V. Ferapontov and K. R. Khusnutdinova, Double waves in multi- dimensional sys- 
tems of hydrodynamic type: the necessary condition for integrability, preprint, arXiv: 
nlin. SI/0412064 

14. A. S. Fokas and L. Y. Sung, On the solvability of the N-wave, Davey-Stewartson and 
Kadomtsev-Petviashvili equations, Inverse Problems 8 (1992), 673-708. 

15. T. Grava and Fei-Ran Tian, The generation, propagation, and extinction of multiphases in 
the KdV zero-dispersion limit, Comm. Pure Appl. Math. 55 (2002), no. 12, 1569-1639. 

16. T. Grava and C. Klein, Numerical solution of the small dispersion limit of Korteweg de 
Vries and Whitham equations, preprint, arXiv: math-ph/0511011 

17. R. H. Hardin and F. D. Tappert, Applications of the split-step Fourier method to the numer- 
ical solution of nonlinear and variable coefficient wave equations, SIAM review, Chronicle 
15 (1973) 423. 

18. T. Y. Hou, J. S. Lowengrub, and M. J. Shelley, Removing the Stiffness from Interfacial 
Flows with Surface Tension, J. Comp. Phys. 114 (1994) 312-338. 



KADOMTSEV-PETVIASHVILI EQUATION 



37 




Figure 35. Solution of the KP-I equation obtained from the ini- 
tial data l|5.1(l with l|5.10jl and v = 10 at time t = 0.4. 



19. G. Huang, V. A. Makarov, and M. G. Velarde, Two-dimensional solitons in Bose-Einstein 
condensates with a disk-shaped trap, Phys. Rev. A 67 (2003), 23604-23616. 

20. E. Infeld, A. Senatorski, and A. A. Skorupski, Numerical simulations of Kadomtsev- 
Petviashvili soliton interactions, Phys. Rev. E 51 (1995), 3183-3191. 

21. E. Infeld, A. A. Skorupski, and G. Rowlands, Instabilities and oscillations of one- and two- 
dimensional Kadomtsev- Petviashvili waves and solitons. II. Linear to nonlinear analysis. 
R. Soc. Lond. Proc. Ser. A (Math. Phys. Eng. Sci.) 458 (2002), no. 2021, 1231-1244. 

22. R. S. Johnson, The classical problem of water waves: a reservoir of integrable and nearly 
integrable equations, J. Nonl. Math. Phys. 10 (2003) 72-92. 

23. C. A. Jones and P. H. Roberts, Motions in a Bose condensate, IV: Axisymmetric solitary 
waves, J. Phys. A Math. Gen. 15 (1982), 2599-2619. 

24. B. B. Kadomtsev and V. I. Petviashvili, On the stability of solitary waves in weakly dis- 
persing media, Sov. Phys. Dokl. 15 (1970), 539-541. 

25. D. Kaya and E. M. Salah, Numerical soliton-like solutions of the potential Kadomtsev- 
Petviashvili equation by the decomposition method, Phys. Lett. A 320 (2003), no. 2-3, 
192-199. 

26. Y. Kodama, A method for solving the dispersionless KP hierarchy and its exact solutions, 
Phys. Lett. A 129 (1988), no 4, 223-226. 

27. Y. Kodama and J. Gibbons, A method for solving the dispersionless KP hierarchy and its 
exact solutions. II, Phys. Lett. A 135 (1989), no. 3, 167-170. 

28. I. M. Krichever, The averaging method for two-dimensional integrable equations, Funkt- 
sional. Anal, i Prilozhen. 22 (1988), no. 3, 37-52; translation in Funct. Anal. Appl. 22 
(1989), no. 3, 200-213. 

29. P. D. Lax and C. D. Levermore, The small dispersion limit of the Korteweg de Vries 
equation 1,11,111, Comm. Pure Appl. Math. 36 (1983), 253-290, 571-593, 809-830. 



C. KLEIN, C. SPARBER, AND P. MARKOWICH 



v = 0.01 V = 0.1 




-6 -4 -2 2 4 -6 -4 -2 2 4 




-6 -4 -2 2 4 -6 -4 -2 2 4 



x x 



Figure 36. Solution of the KP-II equation obtained from the ini- 
tial data 1)5. 1(1 with H5.10JI for several values of v at time t = 0.4. 



30. C. D. Levermore, The hyperbolic nature of the zero dispersion KdV limit, Comm. Part. 
Diff. Equ. 13 (1988), no. 4, 495-514. 

31. D. W. McLaughlin and J. A. Strain, Computing the weak limit of KdV, Comm. Pure Appl. 
Math. 47 (1994), no. 10, 1319-1364. 

32. L. Molinet, J. C. Saut, and N. Tzvetkov, Well-posedness and ill-posedness results for the 
Kadomtsev-Petviashvili-I equation Duke Math. J. 115 (2002), no. 2, 353-384. 

33. S. Novikov, S. V. Manakov, L. P. Pitaevskii, and V. E. Zakharov, Theory of Solitons: The 
Inverse Scattering Method Springer Monographs in Contemporary Mathematics, Springer 
1984. 

34. A. Senatorski and E. Infeld, Simulations of Two- Dimensional Kadomtsev-Petviashvili Soli- 
ton Dynamics in Three- Dimensional Space, Phys. Rev. Lett. 77 (1996), 2855-2858. 

35. G. Schneider, Approximation of the Korteweg-de Vries equation by the nonlinear Schrdinger 
equation, J. Diff. Equ. 147 (1998), no. 2, 333-354. 

36. C. Sulem and P. L. Sulem, The nonlinear Schrddinger equation, Applied Math. Sciences 
139, Springer (1999). 

37. H. Takaoka, Global well-posedness for the Kadomtsev-Petviashvili II equation, Discrete 
Contin. Dynam. Systems 6 (2000), 483-499. 

38. H. Takaoka and N. Tzvetkov, On the local regularity of the Kadomtsev- Petviashvili-II 
equation, Internat. Math. Res. Notices 2001 (2001), no. 2, 77-114. 

39. F. R. Tian, Oscillations of the zero dispersion limit of the Korteweg de Vries equations, 
Comm. Pure App. Math. 46 (1993), 1093-1129. 

40. L. N. Trefethen, Spectral Methods in MATLAB, SIAM, Philadelphia 2000. 

41. http:/ /www. comlab.ox.ac.uk/oucl/work/nick. trefethen 

42. S. Turitsyn and G. Falkovitch, Stability of magneto- elastic solitons and self-focusing of 
sound in antiferromagnets, Soviet Phys. JETP 62 (1985), 146-152. 



KADOMTSEV-PETVIASHVILI EQUATION 



39 



43. N. Tzvetkov, Bilinear estimates related to the KP equations, Journees Equations aux 
Derivees Part. exp. no. XIX, Univ. Nantes, Nantes, 2000. 

44. X. P. Wang, M. Ablowitz, and H. Segur, Wave collapse and instability of solitary waves of 
a generalized Kadomtsev- Petviashvili equation, Physica D 78 (1994), 241-265. 

45. A. Wazwaz, A computational approach to soliton solutions of the Kadomtsev- Petviashvili 
equation, Appl. Math. Comput. 123 (2001), no. 2, 205-217. 

46. V. E. Zakharov and A. E. Kuznetsov, Multi-scale expansion in the theory of systems inte- 
grable by the inverse scattering transform, Physica D 18 (1986), 455-463. 

Max Planck Institute for Mathematics in the Sciences 
E-mail address: kleinSmis.mpg.de 

Wolfgang Pauli Institute Vienna & Faculty of Mathematics, Vienna University, Nord- 
bergstrasse 15, a-1090 vienna, austria 

E-mail address: christof . sparberSunivie . ac . at 

Faculty of Mathematics, Vienna University, Nordbergstrasse 15, A-1090 Vienna, Aus- 
tria 

E-mail address: peter .markowichOuni vie . ac . at 



